function [k] = ShapeK(mth, NmF2, hmF2, B2bot, Az_R)
%
% Shape Parameter
%
%DESCRIPTION:
%This function computes the Shape Parameter (k).
%
%PROTOTYPE:
% [k] = ShapeK(mth, NmF2, hmF2, B2bot, Az_R)
%
%--------------------------------------------------------------------------
% INPUTS:
%   mth        [1x1]       Month                     [-]
%   NmF2       [1x1]       F1-Layer Max. Density     [1e11 m-3]
%   hmF2       [1x1]       F2-Layer Max. Den. Height [km]
%   B2bot      [1x1]       F2-Layer Bottom Thickness [km]
%   Az_R       [1x1]       Eff. Ionis. Level         [-]
%--------------------------------------------------------------------------
% OUTPUTS:
%   k          [1x1]       Shape Parameter           [-]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% 09/02/2022 - Limit on k value added following [2] (pag.14)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Auxiliary Parameter ka
if     (mth == 4) || (mth == 5) || (mth == 6) || (mth == 7)  || (mth == 8)  || (mth == 9)
    ka = 6.705 - 0.014*Az_R - 0.008*hmF2;
elseif (mth == 1) || (mth == 2) || (mth == 3) || (mth == 10) || (mth == 11) || (mth == 12) 
    ka = -7.77 + 0.097*(hmF2/B2bot)^2 + 0.153*NmF2;
end

%Auxiliary Parameter kb
kb = (ka*exp(ka - 2) + 2) / (1 + exp(ka - 2));

%Shape Parameter
k = (8*exp(kb - 8) + kb) / (1 + exp(kb - 8));

% %Alternative Equation [2]
% k = 3.22 - 0.0538*foF2 - 0.00664*hmF2 + 0.113 * hmF2/B2bot + 0.00257*Az_R;

if k < 1 
    k = 1;
end

end